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Abstract:  A  new  approach  to  the  construction  of  a  pos¬ 
teriori  error  estimates  for  finite  element  solutions  of  multi¬ 
parameter  nonlinear  problems  is  presented.  The  esti¬ 
mates  are  derived  from  local,  element-by-element  solu¬ 
tions  of  linearizations  of  the  problems;  they  turn  out  to  be 
very  effective,  computationally  rather  inexpensive,  and 
insensitive  to  the  choice  of  the  local  coordinate  system 
on  the  solution  manifold,  v 


1.  Introduction 


Frequently,  in  practical  computations  in  engineering  and  science,  the  aim  is 
to  obtain  results  which  are  sufficiently  accurate  and  reliable  to  allow  for  a 
decision  about  the  physical  system  under  study.  A  posteriori  error  estimates 
play  a  very  important  role  in  achieving  this  aim.  Such  estimates  are  needed 
not  only  for  judging  the  reliability  of  the  computed  results  but  also  for 
controlling  adaptive  processes  designed  to  achieve  desired  error  tolerances 
at  minimal  cost  or  best  possible  solutions  within  allowable  cost  ranges,  flcf2  ) 

Many  structural  mechanics  problems  are  modelled  by  parameter  depen¬ 
dent  nonlinear  equations.  The  parameters  may  characterize,  for  instance , 
load  points  and  load  directions,  material  properties,  geometrical  data,  etc. 

In  general,  the  set  of  all  solutions  forms  a  differentiable  manifold  in  the  space 
of  the  state  variables  and  parameters.  This  is  often  called  the  equilibrium 
surface  of  the  structure  and  its  form  and  characteristic  features  can  provide 
considerable  insight  into  the  behavior  and  stability  properties  of  the  struc¬ 
ture.  Not  surprisingly,  the  quality  assessment  and  control  of  approximate 
solutions  of  such  parameterized  nonlinear  problems  is  much  more  difficult 
and  expensive  than  that  of  linear  problems.  In  particular,  the  parameter- 
dependence  causes  the  discretization  error  to  become  a  local  concept  which 
depends  on  the  choice  of  the  local  coordinate  system  on  the  equilibrium 
manifold,  [17]. 

In  this  paper,  we  present  a  new  approach  for  computing  a  posteriori 
error  estimates  of  finite  element  solutions  of  nonlinear  equations  involving 

*  This  work  was  supported  in  part  by  ONR-grant  N-00014-90-J-1025  and 
NSF-grant  CCR-8907654. 
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several  real  parameters.  The  approach  is  insensitive  to  the  choice  of  the 
local  coordinate  system  on  the  solution  manifold.  Moreover,  it  has  already 
shown  itself  to  be  highly  effective  and  computationally  inexpensive. 

Many  error  estimations  for  nonlinear  problems  work  with  linearizations 
of  the  given  equations  (see  e.g.  [16]  [19])  and  apply  known  a  posteriori 
error  estimators  for  linear  problems  to  approximate  the  solution-norms  of 
these  linearizations.  So  far,  these  techniques  have  been  very  limited  in  the 
permissible  choices  of  the  local  coordinate  systems  on  the  solution  manifold. 
Another  approach  consists  in  the  defect-based  estimations  proposed,  for 
instance,  in  [10],  [ll],  [12],  [13].  However,  their  widespread  applications 
appears  to  have  been  prevented  by  a  rather  high  computational  cost. 

Our  approach  here  combines  both  approaches.  It  is  based  on  the  use 
of  linearizations  of  the  original  problem,  but  determines  the  required  norms 
of  their  solution  by  solving  the  linearized  equations  only  locally,  element- 
by-element.  In  essence,  the  defect-correction  is  performed  locally  on  the 
linearized  problem.  As  a  result,  the  computational  cost  is  reduced  to  very 
acceptable  levels  that  are  comparable  to  the  cost  of  typical  a  posteriori  esti¬ 
mates  for  linear  problems.  Moreover,  the  resulting  estimates  vary  smoothly 
on  the  solution  manifold  of  the  given  problem  and  apply  to  large  classes  of 
discretizations  and  error  norms. 

In  Section  2,  we  discuss  briefly  some  properties  of  solution  manifolds  of 
parameterized  nonlinear  problems  and  their  discretizations.  Then  Section 
3  presents  the  general  construction  of  the  error  estimates  and  Section  4 
concerns  algorithmic  aspects.  Finally,  some  results  of  numerical  experiments 
for  a  model  problem  are  given  in  Section  5. 

2.  Solution  Manifolds  and  Their  Discretizations 

As  noted  before,  equilibrium  problems  for  many  physical  systems  are  mod¬ 
elled  by  parameter  dependent  nonlinear  equations 

F(x)  —  0,  x  =  (z,\)  (2.1) 

where  z  represents  a  state  variable  and  A  a  vector  of  parameters.  More 
specifically,  suppose  that  the  nonlinear  mapping  F  satisfies  the  condition 

(Fcon)  s  F:ScX*-*Yisa  Fredholm  mapping  of  class  Cr,  r  > 

1  and  index  p  >  1  from  an  open  subset  S  C  X  of  a  real  Banach 
space  X  into  another  such  space  Y. 

Then  it  is  well-known  that  the  set  of  all  regular  solutions, 

M  =  {*;  x  €  5,  F(x)  =  0,  DF(x)X  =  Y}  (2.2) 

is  either  empty  or  a  p-dimensional  CT -manifold  in  X  without  boundary.  We 
assume  always  that  M  /  0 


For  the  numerical  analysis  of  the  solution  manifold  (2.2),  we  need  a 
computationally  feasible  scheme  for  fixing  local  coordinate  systems  at  a 
given  point  x0  €  M.  For  this  suppose  that  a  splitting 

X  =  T  ©  Wt  dimT  =  p  (2.3) 

of  X  has  been  chosen  such  that 

W  nkerDF{x0)  =  {0}.  (2.4) 

Evidently  (Fcon)  and  (2.4)  together  imply  that  the  partial  derivative 
DwF(xq )  is  an  isomorphism  between  W  and  Y.  Now  the  following  result 
(see  [17])  shows  some  neighborhood  of  the  origin  of  T  to  be  diffeomorphic 
to  a  (relative)  neighborhood  of  x0  in  M: 

Theorem  2.1:  Under  the  condition  (Fcon),  suppose  that  at  x0  €  M,  a 
splitting  (2.3),  (2.4)  has  been  given.  Then  there  exist  an  open  ball  B  = 
J3(0,  r)  C  T  centered  at  0  €  T,  an  open  neighborhood  U  C  X  of  x q,  and  a 
unique  Cr -function  r\ :  B  — +  W  such  that  rj( 0)  =  0  and  the  local  coordinate 
mapping 


$  :  B  C  T  t-+  X,  $(a:)  =  x(t)  =  x0  +  t  +  *?(*)»  Vt  €  B  (2.5) 


is  a  CT -diffeomorphism  from  B  onto  M  CiU. 


A  point  x0  €  M  is  called  a  foldpoint  with  respect  to  the  splitting 
(2.3)  if  the  condition  (2.4)  is  violated;  that  is,  if  (2.3)  does  not  induce  a 
local  coordinate  system.  Note  that  the  finite-dimensional  subspace  T  = 
kerDF(xo)  always  admits  a  splitting  (2.3)  for  which  (2.4)  holds. 

In  applications  the  equation  (2.1)  usually  represents  some  boundary 
value  problem  and,  as  indicated  in  (2.1),  we  have  a  natural  parameter  split¬ 
ting  X  =  Z  ®  A,  dimA  =  p  of  X  into  a  state  space  Z  and  p-dimensional 
parameter  space  A.  Evidently,  this  natural  parameter  splitting  may  be  con¬ 
sidered  for  the  definition  of  a  local  coordinate  system.  We  will  assume  that 
there  is  at  least  one  point  x  €  M  where  this  is  possible;  that  is,  where 
Z  fl  ker DF(x)  =  {0}.  Then,  as  noted  earlier,  Z  and  Y  are  isomorphic,  and 
accordingly  there  exists  an  operator  Q  €  L(X,  Y)  such  that  kerQ  =  A  and 
the  restriction  Qq  —  Q\z  €  L{Z,Y)  is  an  isomorphism. 

For  the  computation,  it  is  necessary  to  replace  (2.1)  by  some  finite  di¬ 
mensional  approximating  equation.  Since  the  parameter  space  A  is  already 
finite  dimensional,  only  the  state  space  Z  has  to  be  discretized.  We  follow 
the  approach  in  [17]  in  defining  a  discretization  of  (2.1). 

Let  {11^}  be  a  family  of  finite-rank  projections  IU  €  L(Z)  with  range 
spaces  Zh  =  II \Z,  indexed  by  a  positive  real  number  h  >  0.  Correspond¬ 
ingly  we  introduce  the  extensions  fifc  6  L(X)  defined  by  EU*  =  II 4-  A 
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for  x  =  z  -|-  A  6  X  which  have  the  ranges  Xh  =  Z&  0  A.  With  the  earlier 
isomorphism  Qo  between  Z  and  Y  define  now  the  projections  Pk  €  L(Y), 
Ph  =  Qon^(Qo)-1  with  the  ranges  Yh  =  PhY  =  QoZh-  Then  the  desired 
approximate  equations  are  given  by 

Fh{xh)  =  0,  Xh  =  zh  +  A  €  Xh  (2.6) 

where 

Ph  :  Sh  C  Xh  ►-»  Yh,  Fh(x)  —  PhF(x),  Vx  e  Sh  =  S  Xh-  (2.7) 
Hence,  the  corresponding  discrete  solution  manifold  is 

Mh  =  {x;x  e  Sh,  Fh{x)  =  0,  DFk(x)Xh  =  Ffc}.  (2.8) 

For  the  convergence  theory  additional  assumptions  about  these  dis¬ 
cretizations  are  needed  the  first  of  which  will  be  the  following  consistency 
condition: 

lim  PhV  =  y,  Vy  €  Y.  (2.9) 

Clearly,  any  comparison  of  the  solution  manifolds  M  and  Mh  must 
be  done  locally.  Conceptually,  we  expect  that  for  sufficiently  small  h  the 
local  coordinate  system  at  a  point  x0  €  M  established  by  Theorem  2.1 
also  constitutes  a  local  coordinate  system  at  the  corresponding  approximate 
solution  Xh  €  Mh.  For  this  to  be  correct  we  require  some  stability  condition. 

For  the  analysis  of  the  approximate  problems,  the  discrete  operator 
(2.7)  can  be  extended  to  all  of  S  C  X  as  follows: 

Fh  :  S  Y,  Fh(x)  =  (Jy  —  Ph)Qx  +  PhF(x),  x  €  S.  (2.10) 

Then  the  mentioned  stability  condition  at  xo  €  JVf  assumes  the  form: 

||Z?Ffc(xoMly  ^  7IMIwm  Vw  6  W,  and  sufficiently  small  h  >  0,  (2.11) 

where  7  >  0  is  a  positive  constant  independent  of  w  and  h.  As  shown  in  [17] 
this  condition  is  relatively  easy  to  verify  in  many  practical  situations.  The 
existence  and  convergence  of  solutions  of  the  approximate  problems  (2.6) 
and  the  a  priori  error  estimates  between  M  and  Mh  can  now  be  summarized 
in  the  form  of  the  following  theorem  which  was  proved  in  [17]. 

Theorem  2.2:  (i)  Under  the  condition  (Fcon)  consider  the  discretized 
problem  (2.6)  and  assume  that  the  consistency  condition  (2.9)  holds.  At 
a  given  point  xo  €  M,  suppose  that  the  splitting  (2.3),  (2.4)  induces  the 
local  Cr -coordinate-mapping  (2.5)  and  that  the  stability  condition  (2.11)  is 
satisfied.  Then  for  all  sufficiently  small  h  >  0  there  exists  a  unique  point 
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Xh  €  Mh  and  we  have  lim^o  Xh  =  xq.  Moreover,  the  splitting  X h.  =  T@Wh 
with  Wh  =  W  H  Xh  defines  at  Xh  6  Mh  a  local  Cr -coordinate-mapping 
$h  :  B  C  T  •— ►  Xh-  (ii)  If,  in  addition,  DF  is  Lipschitz  continuous  on 
bounded  sets,  then  there  exists  a  closed  ball  B(0,r0)  C  B,  such  that  for  all 
sufficiently  small  h  >  0  the  estimate 

ll*o(0  -  *fc(t)||*  <  C\\(IY  -  Ph)Qx0(t)\\y,  Vt  €  5(0,  r0), 

holds  with  a  constant  C  that  is  independent  of  h. 

3.  A  Posteriori  Error  Estimates 

For  finite  element  discretizations  of  various  linear  problems  the  theory  and 
application  of  reliable  a  posteriori  error  estimates  has  advanced  considerably 
in  recent  years,  we  refer  only  to  [3]  and  the  two  proceedings  [l],  [5]  where 
many  other  references  are  given.  These  results  are  now  also  being  extended 
to  nonlinear  problems,  see  e.g.  [2],  [4],  [16],  [8],  [19]. 

In  [2]  and  [4]  it  has  been  shown  that  effective  a  posteriori  error  estimates 
and  adaptive  procedures  can  be  constructed  and  incorporated  into  a  general 
continuation  process  for  tracing  paths  on  the  equilibrium  surface.  But  the 
applicability  of  these  results  is  somewhat  restricted  due  to  the  relatively  high 
computational  cost  of  the  estimates.  In  [16]  the  estimates  use  a  linearized 
form  of  the  problem  and  hence  can  be  computed  about  as  rapidly  as  in  the 
case  of  linear  problems.  However,  these  estimates  were  based  on  a  fixed  local 
coordinate  system,  such  as  that  induced  by  the  natural  coordinate  splitting 
X  —  Z  ®  A.  Hence,  by  definition  they  are  not  valid  at  any  foldpoint  of  that 
local  coordinate  system,  and  in  fact,  as  already  simple  examples  indicate, 
the  estimates  may  become  unduly  large  near  any  such  point. 

As  noted  earlier,  our  new  construction  of  a  posteriori  error  estimates  is 
also  utilizes  local  linearizations  of  the  given  nonlinear  mapping.  For  the  re¬ 
lation  between  the  solutions  of  the  nonlinear  and  these  linearized  problems, 
we  present  first  a  simple  result  based  on  tools  from  the  theory  of  Newton’s 
method  (see  e.g.[l4],  [15]) 

Theorem  3.1:  Let  X,  X  be  real  Banach  spaces  and  G  :  S  •— ►  X  a  Cr-map, 
r  >  1,  on  the  open  subset  S  C  X  such  that  DG  is  Lipschitz  continuous  on 
bounded  subsets.  Consider  any  xo  €  S  where  G(x o)  =  0  and  DG(x0)  € 
Isom(X,X).  Then  there  exists  a  closed  ball  B  =  B(x o,r)  C  5,  r  >  0,  such 
that  for  any  x  6  B  the  linearized  problem 

G(x)  +  DG(x)w  =  0  (3.1) 

has  a  unique  solution  w  —  w(x)  €  X  and 

||*o  “  *11  =  (1  +  0(||x  -  xo||))||ts(x)||,  as  x  — +  x0,  x  €  B.  (3.2) 
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Proof.  There  exists  r0  >  0  such  that  B(x0,r0)  C  5  and  DG(x)  is  an 
isomorphism  for  each  x  in  this  ball.  Let  7  be  the  Lipschitz  constant  of  DG 
on  the  ball  and  with  (3  =  ||DG(xo)~1||  and  a  =  /?7  set  r  =  min(r0,  l/(2a)). 
Then  for  all  x  £  B  =  B(xo,r)  the  standard  perturbation  lemma  ensures 
that  DG(x)  €  Isom(X,  X)  and 

l|BC(*r1||  <  t - rp - -jr  <  2/3,  Vi  €  B. 

1  -  a|jx  -  i0|i 

Hence,  w(x)  =  —  DG(x)~1G(x)  is  uniquely  defined  for  x  £  B  and  with 
Xx  =  x  +  w(x)  we  have 

||*i  ~  *o||  <  ||£G(x)-1||||G(®o)  -  G(x)  -  DG(x)(xq  -  x)\\  <  a\\x  -  x0||2, 


Thus 


IM*)ll  <  ll*  -  *o||  +  ||*o  —  *i||  <  (1  +  a||x  -  x0||)||x  -  x0|| 
together  with 

II*  -  *o||  <  ||*  -  *i||  +  ||*i  -  *o||  <  ||tt>(*)||  +  a||*  ~  *o||2 


shows  that 


1  +  a  x  -  xo 


m*)ii  <  11*  -  *oii  < 


a  x  -x0 


M*) 


which  proves  (3.2). 

Now  suppose  that  we  are  in  the  setting  of  both  parts  of  Theorem  2.2.  In 
other  words,  consider  the  problem  (2.1)  where  the  mapping  F  satisfies  the 
condition  (Fcon)  and  DF  is  Lipschitz  continous  on  bounded  sets.  Suppose 
that  a  discretization  has  been  introduced  for  which  the  consistency  condition 
(2.9)  holds.  At  the  currently  given  point  xo  €  M  we  choose  a  splitting 
X  =  T  ®  W,  dimT  =  p  for  which  s0  is  not  a  foldpoint;  that  is,  for  which 
(2.4)  is  satisfied.  Let  r  €  L(X)  denote  the  natural  projection  of  X  onto  T 
along  W.  For  the  mapping 


G:ScX~YxT,  G(x)  =  (F(x),r(x  -  x0)),  Vx  6  5,  (3.3) 


DG(x o)v  =  0  is  equivalent  with  v  £  kerDF(x 0)  n  W  and  thus  v  =  0  and 
DG(x0)  is  injective.  By  (Fcon)  the  mapping  is  also  surjective,  and  hence 
we  find  that  DG(x 0)  6  Isom(X,  Y  x  T)  which  means  that  Theorem  3.1  is 
applicable  to  G. 

Since  xo  €  M  represents  here  the  desired  solution  of  the  original  equa¬ 
tions  (2.1),  we  assume  accordingly  that  x  is  the  corresponding  solution  xj,  of 
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the  discretized  equations  (2.6)  for  which,  by  definition  xo  —  Xh  €  W.  Thus  by 
Theorem  3.1,  the  problem  of  estimating  the  error  xo  —  Xh  is  asymptotically 
equivalent  to  solving  the  constrained  linearized  problem 

F(xh )  +  DF(xh)w  =  0,  t(w)  =  0.  (3.4) 

Of  course,  (3.4)  is  still  an  infinite  dimensional  problem  which  can  only 
be  solved  approximately.  For  this  we  use  the  same  discretization  as  in  the 
computation  of  Xh •  In  other  words,  we  approximate  (3.4)  by  the  constrained 
discretized  problem 

Fh(xh)  +  DFk(xh)wh  =  0,  r(flhwh)  =  0,  (3.5) 

where  the  second  equation  evidently  requires  that  Wh  £  Wh  with  Wh  = 
Xh  H  W.  Because  of  Fh(xh)  =  0,  the  solution  of  (3.5)  is  Wh  =  0.  Hence 
the  norm  ||tr||  of  the  solution  w  of  (3.4)  is  exactly  the  error  arising  in  the 
discretization  (3.5)  of  (3.4).  In  other  words,  by  applying  one  of  the  known, 
linear  a  posteriori  error  estimates  to  (3.5)  we  obtain  asymptotically  valid 
estimates  of  the  discretization  error  j|tu||  and  therefore,  by  Theorem  3.1, 
also  of  the  error  ||xo  —  Xh\\. 

This  was  the  approach  taken  in  [16]  and  [17].  However,  note  that  (3.4) 
and  (3.5)  incorporate  constraints  involving  the  complement  W  of  the  local 
coordinate  space  T.  The  available  a  posteriori  error  estimates  cannot  handle 
arbitrary  choices  of  such  spaces  W\  in  fact,  many  of  them  work  only  when 
W  is  the  given  state  space  Z.  Clearly,  it  is  important  to  allow  for  a  flexible 
choice  of  the  local  coordinate  system  during  the  computation  of  points  on 
the  solution  manifold  of  (2.1).  In  other  words,  this  approach  is  severely 
limitated  by  that  requirement. 

Since  the  discretized  problems  (3.5)  only  produce  the  trivial  approxima¬ 
tion  Wh  =  0  of  the  solution  w  of  the  infinite- dimensional  linearized  problem 
(3.4),  the  question  arises  how  to  obtain  better  approximations  of  w.  Evi¬ 
dently,  the  difficulty  arises  when  we  construct  (3.5)  by  means  of  the  same 
discretization  used  in  the  computation  of  xj,.  In  other  words,  we  should 
work  with  more  accurate  discretizations  of  (3.4).  This  is  the  basic  concept 
of  the  defect  correction  principle.  For  matrix  equations  that  principle  is 
also  called  the  method  of  iterative  improvements  and  was  first  described 
by  Wilkinson  [21].  Extensions  to  ordinary  differential  equations  were  pro¬ 
posed  in  [20]  and  since  then  the  principle  has  been  used  extensively  in  many 
settings,  see  e.g.,  [6],  [7],  [9],  [10],  [11],  [12]  ,  [13]. 

We  shall  not  use  the  full  iterative  form  of  the  defect  correction  principle, 
but  work  instead  only  with  one  specific  improved  discretization  of  (3.4): 

Fh(xh)  +  DF\(xh,)w%  =  0,  wi  €  W*,  (3.6) 
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where  Wh  C  C  W.  In  other  words,  the  system  (3.6)  has  more  degrees 
of  freedom  than  used  in  the  computation  of  x^. 

For  ease  of  notation  we  assume  that  the  spaces  increase  monoton- 
ically  when  h  decreases.  Then  the  convergence  of  the  solutions  wj,  of  (3.6) 
to  that  of  (3.4)  for  h  — *  0  may  be  formulated  in  the  form  of  the  following 
result: 

Theorem  3.2:  Suppose  that  all  conditions  of  Theorem  2.2  hold.  Then 
there  is  some  sufficiently  small  e  such  that  for  allO<h<h<e,  the 
solutions  w  and  of  (3.4)  and  (3.6),  respectively,  exist  and  that 

lim  wi  =  w.  (3.7) 

7L—0 

Proof.  Since  W  and  Y  are  isomorphic,  we  may  choose  some  A  £  Isom(F,  W). 
Then,  by  construction  of  the  local  coordinate  system,  DF(x)A  is  an  iso¬ 
morphism  in  some  ball  Bp  =  B(xo,p)  C  S.  Moreover,  by  the  Lipschitz 
continuity  of  DF  on  bounded  sets,  we  may  choose  p  >  0  small  enough  such 
that  the  stability  condition  (2.11)  holds  for  all  x  £  Bp  and  0  <  h  <  e  with 
a  suitably  small  e  >  0;  in  other  words,  that 

\\DPh(*)Ay\\Y  >  7o||y||y,  Vy  €  Y  (3.8) 

with  some  70  >  0  that  is  independent  of  x  and  h.  Finally,  by  shrinking, 
if  necessary,  p  and  e  further,  we  may  assume  that  for  all  0  <  h  <  e  the 
approximate  solution  s e*.  €  Mh  of  Theorem  2.2  exists  and  that  xh  €  Bp. 

Clearly  (3.8)  implies  that  DFh(x)A  6  Isom(F)  whence  for  any  yi  £  Y 
there  is  a  unique  2/2  €  Y  for  which  DFh(x)Ay2  =  yi-  For  yi  £  Y/,  it  follows 
that 

QAy2  =  PhQAy2  -  PhDF(x)Ay2  +  yi  6  Yh 

and  therefore  Ay2  €  Xh  and  ( lY—Pk)QAy2  =  0;  that  is,  PhDF(x)Ay2  —  yj. 
Hence  we  obtain  from  (3.8)  that 

\\DFh(x)Ay\\y  >  7o||y||y>  Vyen. 

In  other  words,  {DFh(x)A)~r  €  L(Yh)  exists  and  is  uniformly  bounded  by 
6  =  I/70  for  all  x  €  Bp  and  0  <  h  <  e. 

Now  let  h  with  0  <  h  <  t  be  fixed  and  choose  any  h  such  that  0  <  h  <  h. 
By  construction  of  Bp  we  know  that  DF(x\)A  is  an  isomorphism  and  hence 
that  (3.4)  has  a  unique  solution  w  =  Ay  €  W.  Moreover,  the  inverse 
(Z?F^(xfc)i4)-1  €  £(T&)  exists  and  thus  also  (3.6)  has  a  unique  solution 
tt>£  =  Ayi,y  j  €1*.  By  subtraction  we  find  from  (3.4)  and  (3.6)  that 

IK  -  HI  <  *(||(ly  "  Pk)DF(xh)w\ |  +  ||(/y  -  Pk)F(xk) II) 
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which  by  the  consistency  condition  (2.9)  implies  that  (3.7)  holds. 

4.  Algorithm 

As  shown  in  the  previous  section,  for  the  a  posteriori  estimation  of  the 
discretization  error  of  the  computed  point  Xh  G  Mh  we  should  solve  the 
more  accurate  equation  (3.6).  Obviously,  if  X i  =  Xh  then  the  solution 
is  Wk  =  0;  hence  for  any  reasonable  approximation  of  the  discretization 
error  w,  we  require  X-h  to  be  a  sufficiently  larger  subspace  than  Xh-  In 
practice,  the  dimension  of  Xh  is  already  expected  to  be  large  and  thus  a  full 
solution  of  the  equation  (3.6)  will  be  more  expensive  them  the  computation 
of  x h  itself.  This  constitutes  a  severe  restriction  for  the  applicability  of  the 
defect-correction  approach. 

Without  any  further  assumptions  about  the  discretization  process  for 
(2.1)  there  appears  to  be  little  chance  of  resolving  this  problem.  However, 
if  we  suppose  that  (2.1)  is  a  differential  equation  to  which  a  finite  element 
approximation  can  be  applied,  then  it  turns  out  that  an  estimate  of  ||u;^|| 
can  be  obtained  by  solving  (3.6)  locally  on  each  element.  A  detailed  theory 
of  this  approximation  process  would  exceed  the  space  limitations  for  this 
paper  and  will  be  given  elsewhere.  Instead,  we  shall  sketch  here  only  the 
general  computational  procedure. 

We  consider  a  nonlinear  problem  in  a  generic  weak  form  requiring  the 
determination  of  z  £  Z  such  that  for  given  A  6  A 

b(z,  A,u)  =  g(\,v),  VveV  (4.1) 

where  Z  and  V  are  real  Hilbert  spaces  and  A  is  a  p-dimensional  inner- 
product  parameter-space.  In  (4.1),  b  stands  for  a  form  on  Z  x  A  x  V  which 
is  nonlinear  in  the  first  two  variables  but  linear  in  the  third  one,  and  g  is 
a  functional  on  A  x  V  which  is  linear  in  v  but  may  be  nonlinear  in  the 
parameter-variable.  We  introduce  also  the  Hilbert  space  X  =  Z  x  A  with 
the  usual  inner  product  for  product  spaces  denoted  by  (.,.). 

The  discretizations  are  now  specified  by  the  choice  of  finite  dimensional 
linear  subspaces  Zh  C  Z  and  C  V  and  the  corresponding  orthogonal 
projections  II ^  €  L(Z)  onto  Zh-  In  other  words,  the  resulting  (Galerkin- 
type)  discretization  (2.6)  becomes  the  problem 

determine  Zh  €  Zh  such  that 

b(zh,\,Vh)  =  g{ A,vfc),  vh  G  Vh.  (4.2) 

If  z^,  j  =  1, 2, . . . , mh  and  vJh,  j  =  1,2,..., m/,  denote  bases  of  Zh  and  Vh, 
respectively,  then  (4.2)  assumes  the  usual  form 

b(£,Cizh’X’vh)  =9(*>Vh),  k  =  l,2,...,mh  (4.3) 

;=i 
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of  a  system  of  m,h  nonlinear  equations  in  -f-  p  variables. 

At  any  solution  x  =  ( z ,  A)  of  (4.1)  consider  now  a  local  coordinate  sys¬ 
tem  defined  by  a  p- dimensional  subspace  T  C  X  together  with  its  orthog¬ 
onal  complement  W  =  T1-  in  X.  Under  appropriate  smoothness  assump¬ 
tions  on  b  and  g  the  linearized  problem  (3.4)  at  a  computed  approximation 
xx  €  Xh  =  Zh.  ©  A  of  x  then  becomes 

determine  u  €  Z,  p  £  A  such  that 

ah(u,p,v)  =  f(v),  Vu  £  V,  (4.4) 

(tk,(u,p))  =  0,  fc  =  l,...,p 

where  t*  EX,  k  =  1,2, . . .  ,p  represents  a  basis  of  T  and 

ah(u,/j.,v)  =  Dzb(zh,  A h,v)u  +  D\b(zh,  \h,v)p  -  Dxg(Xh,v)p  (4.5) 
f(v)  =  -b(zh,  Xh,v)  +g{\h,v) 

As  we  saw,  we  have  to  solve  (4.4)  with  a  discretization  induced  by  some 
larger  subspaces  Z%  =  Zh  ©  Z£  and  14  =  14  ©  V£  where  Zch  and  V£  denote 
here  certain  complements  of  Zx  and  14  in  Zj^  and  V),  ,resp  ectively.  Hence, 
the  resulting  discretization  is 

determine  Uh  £  Zh,  uch  £  Z{,  p  6  A  such  that 


ah(uh,P,vh)  =  f(vh),  €  14,  (4.6a) 

ah(ul,p,Vh)  =  0,  Vvh  £  14,  (4.6b) 

ah(uh  +  ueh,p,vch)=  f(veh),  £  V£,  (4.6c) 

(ik,(uh,p))=  0,  k=l,...,p  (4 .6d) 

(**»«>/*))  =  0,  k=l,...,p  (4.6e) 


But  (4.6a)  and  (4.6d)  constitute  exactly  the  discretization  of  (4.4)  defined 
by  the  restriction  ux  €  Zh  and  Vh  £  14,  for  which  we  found  the  solution  to 
be  zero.  Hence,  it  remains  only  to  solve  the  reduced  problem 

determine  uj  E  Z(,  p  £  A  such  that 
ah(uch,p,vch)=  f(v{),  Vr£6l4c, 
ah(uch,p,vh)  =  0,  £  14, 

fc=l,...,p 

where  (4.7b)  represents  a  boundary  condition  for  uch. 

For  the  practical  application  we  have  to  choose  a  suitable  space  T. 
Since  the  computed  point  Xh  belongs  to  the  manifold  Mh  of  the  discretized 
problem  (4.3)  it  is  natural  to  select  T  as  the  null-space  of  the  Jacobian 


(4.7a) 

(4.71.) 

(4.7c) 
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of(4.2);  that  is,  as  the  subspace  of  Xh  corresponding  to  the  tangent  vector 
of  Mh.  at  Xh-  If  A  is  one-dimensional  and  Mh  is  computed  by  a  standard 
continuation  process  then  a  normalized  basis  vector  of  T  is  usually  available 
at  each  computed  point.  Analogously,  in  the  multi-parameter  case,  if  a 
triangulation  of  Mh  is  being  computed  (see  e.g.  [18])  then  again  orthonormal 
"tangent-bases”  tJh,  j  =  1, . . .  ,p  axe  readily  available  at  the  computed  points 
on  the  manifold. 

Let  ft*  be  a  typical  element  of  the  mesh  on  the  underlying  domain  ft 
of  the  original  problem  (4.1)  corresponding  to  the  finite  element  space  Zh- 
In  order  to  obtain  D  Zh,  we  may  subdivide  ft*  into  sub-elements  or 
increase  the  order  of  the  element  or  both,  in  other  words  we  may  use  either 
a  local  h-version,  p-version,  or  h-p- version  of  the  finite  element  method. 
This  has  the  advantage  of  being  applicable  to  large  classes  of  error  norms 
for  u.  Of  course,  for  the  solution  of  (4.7a, b,c)  it  now  will  become  necessary 
to  "extend”  the  computed  point  Xh  and  its  corresponding  "tangent "-basis 
tJh,  j  =  1, ...  ,p  to  by  suitable  interpolation.  The  norm  of  the  computed 
local  solution  on  ft*  is  called  the  local  error  indicator  on  that  element;  that 
is, 

Vj  =  IKIIqk,  (4.8) 

and  the  error  estimate  for  the  computed  point  Xh  is  the  sum 

*  =  (I>i)1/2.  (4-9) 

i 

over  the  error-indicators  of  all  elements  of  ft. 

Altogether,  for  the  computation  of  the  error  estimate  at  any  computed 
point  of  the  solution  manifold  Mh  of  the  discreterized  problem  we  have  now 
the  following  algorithm: 

1.  Let  Xh  €  Xh  on  Mh  be  the  current  solution  of  (4.2).  Select  the  tangent 
space  Th  =  span{t\, . . .  ,t£}  of  Mh  at  Xh  to  define  the  local  coordinate 
system. 

2.  Loop  over  all  elements  ft*  of  the  domain  ft 

2.1  Subdivide  the  element  ft*  into  sub-elements  corresponding  to  the 

desired  choice  of  a  finite  dimensional  subspace  Zj[  on  ft*. 

2.2  On  the  refined  element  solve  the  local  linearized  finite  element 
problem  (4.7a,b,c)  for  (uch,p)  G  X^  =  Zj!  0  A.  This  involves  the 
interpolation  of  x\  and  t*,  •  •  • ,  t*  from  X3h  into  X^. 

2.3  Compute  the  local  error  indicator  (4.8)  for  the  element. 

3.  Compute  the  desired  a  posteriori  estimate  (4.9)  for  the  current  point 
Xh  €  Mh' 

5.  Numerical  Example 


ll 


As  a  model  problem  for  numerical  experiments  we  consider  the  two-dimen¬ 
sional  nonlinear  boundary  value  problem 


-A2  =  Aez,  z  =  z(x,y),  V(z,  y)  €  ft  =  (0, 1)  x  (0, 1),  (5.1) 

2  =  0  on  dfl. 


A  weak  formulation  is 

I  (zxvx  +  ZyVy  -  \ezv)dxdy  =  0,  Vu  €  Hq(CI).  (5.2) 

J  n 

and  we  assume  that  A  €  R1  which  means  that  the  solutions  of  (5.2)  are 
expected  to  form  a  one-dimensional  manifold  M. 

For  the  discretization  we  use  a  uniform  mesh  of  sixteen  bi-quadratic 
elements  on  ft.  For  the  computation  of  the  one-dimensional  solution  man¬ 
ifold  Mh  of  the  discretized  problem  a  continuation  process  (PITCON)  is 
applied  starting  from  (z°,  A°  )=  (0,0)  and  our  aim  is  to  determine  a  pos¬ 
teriori  estimates  of  the  error  between  M  and  Mh  at  all  computed  solutions 

(^/ii  A/i)  £  Affc. 

In  line  with  (4.4),  the  linearized  problem  at  ( z ,  A)  has  the  form 

J  [wiU*  +  uyvy  —  (Au  •+-  y)ezv\  dx  dy 
n 

=  —  [  (zxvx  +  ZyVy  -  \ezv)dxdy,  Vv  6  (5.3) 

Jn 

As  before  the  point  (2,  A)  stands  here  always  for  one  of  the  computed  points 
(zh.i  Afc)  €  Mh-  We  apply  now  the  algorithm  of  Section  4.  and  approximate 
the  error  ||u||  by  solving  (5.3)  locally  on  a  larger  finite  element  subspace  Z\ 
of  J5T$(ft)  together  with  the  auxiliary  condition 

(**,«)  =  0,  (5.4) 


Here,  as  noted  in  Section  4.,  we  choose  th  as  a  normalized  tangent  vector  on 
Mh  at  (zh,  Afc).  Such  a  tangent  vector  is  available  at  each  step  of  the  con¬ 
tinuation  process  and  hence  the  equation  (5.4)  will  involve  little  additional 
computational  cost. 

For  the  local  solution  each  one  of  the  16  elements  of  ft  were  divided  into 
(k  4-  l)2  biquadratic  sub-elements  with  k  =  1,2, 3, 4.  The  resulting  error 
estimates  are  shown  in  Table  5.1  where  |ju^fc||  denotes  the  computed  error 
norms  for  the  four  cases  of  k.  The  computations  are  very  cost-effective, 
since  each  local  problem  involves  only  a  fixed  number  of  degrees  of  freedom 
depending  on  the  value  of  k.  If  ||u&4||  is  taken  as  the  exact  error  then 
the  effectivity  index  of  the  estimates  is  surprisingly  good  even  for  k  =  1. 
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The  table  also  shows  that,  as  expected,  the  estimated  errors  vary  smoothly 
along  the  solution  path  Mk  and  show  no  sudden  increases  near  the  limit 
point  A  =  6.804524. 

As  mentioned  earlier,  if  the  natural  coordinate  system  induced  by  the 
parameter  space  A  is  chosen  then  we  expect  the  resulting  error  estimates 
to  become  unduly  large  near  the  limit  point.  This  is  indeed  the  case  as  the 
last  column  of  Table  5.2  shows.  At  the  same  time,  it  should  be  noted  that 
the  computational  cost  of  the  two  approaches  are  practically  identical. 
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Table  5.1 


step 

A 

KJi 

■ran 

IK,  II 

mm 

■ 

1 

0.118710 

0.000516 

0.000538 

0.000542 

0.000544 

2 

0.473781 

0.002040 

0.002128 

0.002146 

0.002152 

3 

0.946226 

0.004010 

0.004184 

0.004220 

0.004231 

4 

1.415882 

0.005898 

0.006155 

0.006208 

0.006225 

5 

1.882269 

0.007695 

0.008031 

0.008101 

0.008123 

6 

2.344791 

0.009389 

0.009801 

0.009887 

0.009914 

7 

2.785973 

0.010913 

0.011395 

0.011495 

0.011527 

8 

3.220584 

0.012313 

0.012861 

0.012976 

0.013011 

9 

3.647333 

0.013578 

0.014188 

0.014315 

0.014355 

10 

4.064580 

0.014694 

0.015361 

0.015500 

0.015544 

11 

4.470212 

0.015648 

0.016368 

0.016519 

0.016555 

12 

4.861489 

0.016433 

0.017200 

0.017362 

0.017413 

13 

5.234830 

0.017050 

0.017861 

0.018032 

0.018087 

14 

5.585562 

0.017522 

0.018373 

0.018554 

0.018611 

15 

5.907655 

0.017914 

0.018805 

0.018994 

0.019054 

16 

6.193552 

0.018368 

0.019300 

0.019499 

0.019562 

17 

6.434373 

0.019144 

0.020128 

0.020337 

0.020404 

18 

6.620862 

0.020640 

0.021692 

0.021916 

0.021986 

19 

6.745397 

0.023317 

0.024468 

0.024710 

0.024787 

20 

6.804524 

0.027539 

0.028832 

0.029100 

0.029185 

21 

6.800451 

0.033461 

0.034945 

0.035250 

0.035345 

22 

6.740221 

0.041066 

0.042795 

0.043146 

0.043255 

23 

6.633252 

0.050283 

0.052305 

0.052712 

0.052838 

24 

6.489009 

0.061065 

0.063424 

0.063896 

0.064041 

25 

6.328924 

0.072480 

mrnsmm 

0.075729 

0.075894 

26 

6.144920 

0.085598 

0.088703 

0.089318 

0.089506 

27 

5.941748 

0.104046 

0.104745 

■»n  iimu 

28 

5.723487 

0.117289 

0.121324 

0.122117 

0.122359 

29 

5.493577 

0.136117 

0.140694 

0.141591 

0.141865 

30 

5.254879 

0.157154 

0.162338 

0.163350 

0.163659 

Table  5.2 


step 

A 

KJI 

INI 

15 

5.907655 

0.019054 

0.022905 

16 

6.193552 

0.019562 

0.024439 

17 

6.434373 

0.020404 

0.027792 

18 

6.620862 

0.021986 

0.041992 

19 

6.745397 

0.024787 

0.159595 

20 

6.804524 

0.029185 

0.176171 

21 

6.800451 

0.035345 

0.086353 

22 

6.740221 

0.043255 

0.069906 

23 

6.633252 

0.052838 

0.065171 

24 

6.489009 

0.064041 

0.065099 

25 

6.328924 

0.075894 

0.067672 
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